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Abstract 


This study aims to present an accelerated derivative-free method for solv- 
ing systems of nonlinear equations using a double direction approach. The 
approach approximates the Jacobian using a suitably formed diagonal ma- 
trix by applying the acceleration parameter. Moreover, a norm descent 
line search is employed in the scheme to compute the optimal step length. 
Under the primary conditions, the proposed method’s global convergence 
is proved. Numerical results are recorded in this paper using a set of 
large-scale test problems. Moreover, the new method is successfully used 
to address the problem of Chandrasekhar’s integral equation problem ap- 
pearing in radiative heat transfer. This method outperforms the existing 
Newton and inexact double step length methods. 
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1 Introduction 


Scientists are interested in nonlinear problems because most engineering, bi- 
ology, mathematics, physics, and other science problems are naturally non- 
linear. The standard nonlinear equation system is represented by 


F(x) =0, (1) 


where F : R” — R” is a nonlinear map. The space R” denotes the n- 
dimensional real space, || - || is the Euclidean norm, and Fy = F(x,) is used 
throughout this paper. Further applications of problem (1) can be found 
in chemical equilibrium systems [16] and signal and image processing [25]. 
The Chandrasekhar H-equation that arises in the theory of radioactive heat 
transfer is a nonlinear integral equation that can be discretized into non- 
linear equations [20]. Iterative methods for solving these problems include 
the Newton and quasi-Newton methods [4, 21, 26, 14], Levenberg—Marquardt 
methods [15, 13, 12], matrix-free methods [17, 1, 8], and tensor methods [2]. 
Typically, the iterative formula for solving these methods is given by 


Ce+1 = Xe +argdy, k=O0,1,..., (2) 


where x41 represents a current iterate, x; is the previous iterate, a; is a step 
length, and d; is the search direction can be calculated by solving system of 
linear equations as follows: 


Fy, + Fidp = 0, (3) 


where FY, is the Jacobian matrix of Fy, at x,. One of the most important 
requirements of the line search is to reduce the function values sufficiently 
[9, 11], as shown below: 

I|Fetill < |Fell- (4) 


Irrespective of how appealing the Newton and quasi-Newton approaches 
are, the Jacobian matrix or its approximation can be calculated at each 
iteration, making them unsuitable for solving large-scale problems. Due to 
the drawbacks of these methods, the double direction technique has been 
proposed [6], with the following iterate: 


Chri = ty + apd, + a2 dz, (5) 
where dx and by are search directions, respectively. 


Suppose that f is a merit function defined by 


fle) = SIF@P, (6) 
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Then problem (1) is analogous to the unconstrained optimization problem 
described below: 
min f(z), «2 €R", (7) 


where f : R” > R, and condition (4) is equivalent to 
F(t + andy) < f (zx). (8) 


The iterative method generating the sequence {x,} that satisfies (4) is called 
the norm descent method [9]. If d;, is a descent direction of f at x,, then con- 
dition (8) holds for all a, > 0 small enough. The Newton method (NM) with 
line search is norm descent. Nonetheless, d, might not be a descent direction 
of f at x, for quasi-Newton methods, even if the approximation of the Ja- 
cobian matrix By is positive definite and symmetric. Li and Fukushima [14] 
proposed an approximately norm-descent line search approach. However, the 
proposed method is not norm descent, but they established a global conver- 
gence theorem under the assumption that Jacobian is uniformly nonsingular. 


The concept of the double direction approach was suggested by Duranovié- 
Miliéié [6] by using a multi-step iterative scheme and curve search to generate 
new iterates. However, in [7], another double direction algorithm was also 
presented to minimize nondifferentiable functions. Motivated by the work 
presented in [7], Petrovié and Stanimirovié [19] suggested a double direc- 
tion model for solving unconstrained optimization problems. They used the 
acceleration parameter 7, to approximate the Hessian matrix, that is, 


Vf (xe) © Vel, (9) 


where J is the identity matrix, and the sequence of iterates {a,} is gener- 
ated using (5). The attractive feature of the scheme in [19] is that the two 
directions presented in their work are derivative-free. Therefore, it enables 
their method to solve large-scale problems. However, the literature is infre- 
quent to study derivative-free double direction methods for solving nonlinear 
equations. Based upon the idea presented in [19], Halilu and Waziri used 
the scheme in (5) to propose a method for solving a system of nonlinear 
equations using a double direction approach. They used the acceleration 
parameter 7, > 0 in their work to approximate the Jacobian matrix, that is, 


Fw el, (10) 


where J is an identity matrix and the acceleration parameter is derived as 


Va Yk (11) 


Vk+1 = : 
(an + of yn) yb de 


The method’s global convergent is proved by assuming that the Jacobian of 
F is positive definite and bounded. The double direction scheme is justified 
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by the fact that scheme (5) contains two corrections. If one of the iterative 
corrections fails, then the system will be corrected by the second. 

The implementation of double direction is additionally enhanced by 
Petrovié [18], where the double step length scheme for the unconstrained 
optimization problem is presented as 


Lk+1 = Le + Ande + Bebx, (12) 


where a, and {, are two different step lengths. The numerical results indi- 
cated the approach is quite effective compared to the double direction method 
in [19]. The authors in [8] incorporated the concept in (12) and transformed 
the double step length method for solving (1) to improve the numerical re- 
sults and global convergence properties of the double direction scheme. The 
numerical results exhibited that the method in [8] is more reasonable than 
the method in [10] because it converges faster. Furthermore, the method 
[8] is globally converged using the line search proposed in [14]. Motivated 
by the work in [8], Halilu and Waziri [11] presented an inexact double step 
length method for solving (1). The attractive feature of this method is that 
it has a double step length and a single direction that satisfies the decent 
properties independent of line search. Despite the good convergence prop- 
erties of the method in [10], its numerical performance is defined as weaker. 
Therefore, motivated by this reason, we aim to develop a globally converged 
derivative-free method with a line search to solve a system of nonlinear equa- 
tions without calculating the Jacobian matrix. 


Table 1: Authors’ contribution table 


Author’s Name Derivative-free | Matrix-free | Double Direction | Global Convergence | Application 

Duranovié-Miliéié [6] v v v v 
Halilu and Waziri [10] v v v v 
Duranovic-Milicic [7] v v v v 
Musa, Waziri, and Halilu [17] v v 

Abdullahi, Waziri, and Halilu [1] v v v v 
Petrovié and Stanimirovié [19] v v v v 
Kanzow et al. [12 v v 
Halilu and Waziri [11] V v v v v 
Yuan and Lu [26 v v 

Waziri et al. [23 v v v 
Halilu and Waziri [8] v v v v 
Li and Fukushima [14] v v 
Halilu and Waziri [9] v v v v 
petrovié [18] v v v v 

This article v v v v v 


The research gap between the existing method and this article is described 
in Table 1 above. The table clearly shows that only the proposed method 
is derivative-free, matrix-free, double direction method, globally convergent, 
and can be applied to solve discretized Chandrasekhar’s integral equation 
among the listed articles. 

We now describe how the paper is structured. The proposed method’s 
algorithm will be presented in section 2. Section 3 illustrates the convergence 
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results. Section 4 contains a list of numerical experiments and applications 
of the proposed method to Chandrasekhar’s integral equation, which arises 
in radiative heat transfer. Section 5 concludes the paper. 


2 Main result 


In this section, we present the algorithm of our method. We suggest that the 
dy and by in (5) are defined as follows: 
di, = 1, Fe (13) 


and 
by = —F, (14) 
where y, > 0 is an acceleration parameter. By substituting (13) and (14) 


into (5), we obtain 


Uk+1 = LE (ay + OR VK) Vy Fe. (15) 
The acceleration parameter can be obtained using Taylor’s series expansion 
below: 
Frat ~ Pet F'(w) (rey — £k). (16) 
By multiplying (16) through by 6;, we have 
On Fea. © OnFR + O4.F" (a) (arr + aryn)de, (17) 
where 6, > 0 and w satisfies the conditions w € [ap, e~41] and 
p=tetC(ee+1- ee), OSC <1. (18) 


Taking ¢ = 1 in (18), obtain = p41. 
We like to make Jacobian approximations via 


O,F' (ab) & Yel. (19) 


Using (17) and (19), it is easy to confirm that 


YVe+1Sk = IYK, (20) 
aM 
where s, = (ax +07YK)dk, Ye = Feti — Fy, and 0, = Ts, (see [24]). 
kok 


The proposed acceleration parameter is defined by multiplying y? on both 
sides of (20) 
Il sell? lvl? 


= . 21 
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Our proposed scheme is given by equation (22) below based on (13) and (15): 


Chri = te + (ag + 027%) dg. (22) 


Algorithm 1 Modification of the double direction approach (MDFDD) 
Input: Given x9, Yo = 1, € = 10~°, d; > 0, d2 > 0, and r € (0,1), set k = 0. 
Step 1: Compute Fy. 

Step 2: If || F;|| <, then stop; otherwise, proceed to Step 3. 

Step 3: Calculate search direction d, = —, Fr. 

Step 4: Set rp41 = pe + (ax +07 7%)de, where a, = 7% with az being the 
smallest nonnegative integer a such that 


f (te + (an +0%%e) de) — f (te) < —dbr\lonFel|? — ¢2llande ||? + tef (we). (23) 


Let {7} be a given positive sequence such that 
Co 
bz Th <T< OO. (24) 
k=0 


Step 5: Compute Fy41. 
Il sell? llvell? 


(an + ORK)? (YE die)? 
Step 7: Consider k = k+ 1 and go to Step 2. 


Step 6: Determine 7x41 = 


3 Convergence Analysis 


We present how the proposed Algorithm 2 (MDFDD) converges globally in 
this section. Let us start by defining the level set 


Q = {2|||F(x)|] < ||F(@o)|[F- (25) 
However, we require the following assumptions: 
Assumption 1. However, we state the following assumptions: 
1. There exists «* € R” such that F(a*) = 0. 


2. F is continuously differentiable in some neighborhood say Q of «* con- 
taining 2. 


3. The Jacobian of F' is bounded and positive definite on Q. That is, there 
exist positive constants H > h > 0 such that 


||F’(2)||< H forallze€Q, (26) 
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and 
All|? < d?F'(x)d_ for allx €Q,d€ R”. (27) 


Remark 1. We make the following remark: 
Assumption 1 implies that there exist constants H > h > 0 such that 


Alldl| < |F’(a)d|| < H|d||_ for allxeQ,deR”, (28) 
Alley SF@)-FO)|< Alla—yl| foralla,yeQ. (29) 


Since yJ approximates Fy along sz, the following assumption can be 
made. 


Assumption 2. 7;J is a good approximation to F7, that is, 
(Fa — YL) dell < ellFell, (30) 
where € € (0,1) is a small quantity [26]. 


Lemma 1. Suppose that Assumption 2 holds, and let {x;,} be generated by 
the MDFDD algorithm. Then d;, is a sufficient descent direction for f(x,) 
at rz, that is, 

Vf (ar)* dy <¢l|Fil|?, ¢> 0. (31) 


Proof. From (13), we have 
Vi (ax)? dy = FE Fidx 
= Fe [Fi — vel) dn — Fx] (32) 
= FE (Fe — Yel) de — || Fell?, 
by the Cauchy—Schwarz inequality, we have 


VF (ae)? de < | Fellll(Fe — WL) dell — F(x) IP? 


33 
< -(1- | F(ax)|I?. oe 


This lemma is true for ¢ € (0,1). 
We can conclude from Lemma 1 that the norm function f(a;) is a descent 
along d;,, which means that || Fy+41|| < || F%|| is true. 


Lemma 2. Suppose that Assumption 1 holds, and let {x,} be generated by 
the MDFDD algorithm. Then {x} CQ. 


Proof. From Lemma 1, we have ||Fy+1|| < || £;||. Furthermore, for all k, 


|Fe+all < ||Fell < |]Fe-all < +++ < ||Foll- 


This means that {x,} CQ. 
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Lemma 3 (see [26]). Suppose that Assumption 1 holds, and let {a;,} be 
generated by the MDFDD algorithm. Then there exists a constant m > 0 
such that for all k, 

UE Sk > hllsell?. (34) 


Lemma 4. Suppose that Assumption 1 holds and that {2,} is generated by 
the MDFDD algorithm. Then 


lim |/a,dx|| = lim ||s,|| = 0 (35) 
k-00 k- 00 
and 
k—- 00 


Proof. By (23) for all k > 0 


Pallaxdal|” < dillon Fell? + d2lloxdell? 


(37) 

< || Fall? — || Fetal? + rll Fall?. 

By summing the above inequality, we have 
k k k 
2S > llasdall? < $5 (Fill? — Final?) + 0 nl LAA, 
i=0 i=0 1=0 
k 

= ||Foll? - Feu? + 30 FP, 

= (38) 


IA 


k 
Fo||? + Full? S > 7, 
i=0 


CO 
< ||Foll? + Full? 
i=0 


From the level set and the fact that {7,} satisfies (24), then the series 


S- \|a;d;||? converges. This implies (35). Using the same logic as above, 
i=0 
but this time with ¢,||a,F%||? on the left, we obtain (36). 


Lemma 5. Suppose that Assumption 1 holds, and let {x,} be generated by 
the MDFDD algorithm. Then there exists a constant m; > 0 such that for 
alk > 0, 

I|dx || < B. (39) 


Proof. From (13) and (21), we have 
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(ye _ yeni)? F 
Iyer? [1s%—1112 
|| Fell lse—a1/2lye—all? (40) 


Ilsx—1 11? [yea 11? 
< ||Fol|- 


Ilda] = 


IA 


Choosing B = ||Fo||, we have (39). 


Theorem 1. Suppose that Assumption 1 holds, and let {z,} be generated 
by the MDFDD algorithm. Assume further, for all k > 0, 


|Fie de| 
k 2 r 5) (41) 
I| dx ||? 
where A is some positive constant. Then 
lim || F),|| = 0. (42) 
k-0o 


Proof. From Lemma 5, we have (39). Also, from (35) and the boundedness 
of {||dx||}, we have 


lim ag||dy|? = 0. (43) 
k—-00 
From (41) and (43), we have 
lim |Fod,| = 0. (44) 
k—-00 
Also, from (13), we have 
Fe dy, = — 75 (Fell? (45) 
|| Fell? = || — Fe devel (46) 
< |Fe dello 
Since 
yo = (ye—18k—1)” ay h?||sx—1||* S h? || sx—1|? _ h? 
Myer ilse—all? ~ Mya IP Isai? ~ HP Ilse—all? A?’ 
then 
= h? 
ln | 2 2’ 


Therefore from (46), we have 


H? 
IFIP < IF dl (i) (47) 
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As a result, 
2 T He 
0< |Fil? < IFEaal (Fe) 0. (48) 
Hence, 
lim ||F;|| = 0. (49) 
k—- oo 


4 Numerical experiments 


The first part of this section provides numerical results to demonstrate the 
efficacy of the proposed method by comparing it with existing methods. 


e IDFDD: Algorithm 1 proposed in [10]. 
e IDSL: Algorithm 1 proposed in [11]. 


The MDFDD method is used in the second part to solve the problem of 
Chandrasekhar’s integral equation in radiative heat transfer. The computer 
codes used were written in MATLAB 9.4.0 (R2018a) and ran on a computer 
with a 1.80 GHz CPU processor and 8 GB RAM. 


4.1 Experiment of some nonlinear systems of equations 


In the experiments, we implemented the three algorithms using the same 
line search (23), with ¢, = ¢2 = 107+, r = 0.2, and ™% = wa The 
iteration is set to stop for the three methods if ||F;,|| < 10~° or when the 
number of iterations overreaches 1000, but there is no x, meeting the stopping 
criterion. The numerical effects of the three methods are shown in Tables 3-9, 
where “ITRN,” “CTM(S),” and “IP” represent the total number of iterations, 
CPU time (in seconds), and initial points, respectively. In addition, ||F’&|| 
represents the residual value at the stopping point. The symbol “-” indicates 
failure due to a memory requirement or when some iterations exceed 1000. 
We tested the three methods on the current seven test problems, each with a 
different set of initial points and dimensions (n values). The experiment was 
carried out with the dimensions 100, 1,000, 2,000, 10,000, 50,000, and 100,000 
to demonstrate the comprehensive numerical experiments of the MDFDD, 
IDFDD, and IDSL methods. Table 2 contains the starting points for the test 
problems. 
The experiments made use of the following test problems: 

Problem 1 [8] 


cos( #1422 ) 


n+1 


Fi =2z,-e 
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Table 2: Initial points used in test problems 


INITIAL POINTS (IP) VALUES 


IP1 ees 

IP2 ee a 

IP3 (3,3,053)- 

IP4 Cae 

IP5 O42 ond) 
n\T 

IP6 (4.44... 4") 

IP7 (ee 


@i-1teitei4y 
n+l 


cos (“2=4 *) 
a aT — 
Fy, = %n—e ne , t=2,3,...,n—1. 


Problem 2 [11] 
F,(x) = 2j(1 + 242n-2%n-10n) —-2+ (1-27), i=1,2,...,n. 


Problem 3 [24] 


Problem 4 [10] 
Fy (2) = (aj + 23)21 — 1, 
Fi(x) = (wy + 27 + @741)ei — 1, 


Fi@a=@a +e )e,, 1=2)3,4.0—1. 


Problem 5 [10] 
F(x) = 2x;—sin|x;|, i=1,2,...,n. 


Problem 6 [11] 
F(a) = Ax + by, 


where A = ia Sg. te , and by = (e% —1,...,e% — 1)”. 


Problem 7 [8] 
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F(x) = Ba + bo, 

2-1 

02 -1 
where B = a Soon A , and be = (sinz; —1,...,sinz, —1)?. 
area 
—1 2 


Table 3: Numerical outcomes of MDFDD, IDFDD, and IDSL methods for 
problem 1 


MDFDD IDFDD IDSL 
Dimension IP ~ITRN  CTM(S) F,|| “TERN CTM() F,|| “ITRN CTM(S \| Fe 
100 IP 6 0.009 7.44E-06 60 0.031  8.35E-06 I 0.015  9.86E-06 
IP2 0.020 267E-06 60 0.023 9.48E-06 42 0.013 7.83E-06 
IP3 6 0.013. 9.57E-06 58 0.022 7.93E-06 40 0.022 7.72E-06 
IP 7 0.009 7.55E-06 60 0.020 8.73E-06 42 0.015 7.21E-06 
P57 0.023 7.17E-06 59 0.023 —-8.78E-06 1 0.014 7.88E-06 
IP6 8 0.024 7.36E-06 61 0.018 849E-06 42 0.026 — 9.23E-06 
IP7 10 0.027 668E-06 61 0.019  7.64E-06 42 0.023 8.3E-06 
1,000. +IP 3 0.005 435E-07 96 0.057 859E-06 45 0.031 7.51E-06 
P23 0.005 5.13E-07 96 0.077 9.76E-06 45 0.051 8.52E-06 
IP3 3 0.007 2.09E-07 94 0.056 817E-06 43 0.029 8.41 E-06 
IP 3 0.005 46E-07 96 0.058  8.98E-06 45 0.057 -7.84E-06 
IP5 3 0.009 3.17E-07 95 0.057 88E-06 44 0.035 8.34E-06 
IP6—3 0.007 6.38E-07 97 0.058 8.74E-06 46 0.055 7.03E-06 
IP7 3 0.011 5.66E-07 97 0.092 7.98E-06 45 0.051 9.17E-06 
10,000 IP 3 0.017‘ 1.38E-10 1 0486 9.05E-06 48 0.239 8.14E-06 
IP2 3 0.035 — 1.63E-10 2 0.492 7.81E-06 48 0.206 —-9.24E-06 
IP3 3 0.022 664E-11 109 0.477  8.6E-06 6 0.266 — 9.13E-06 
IP 3 0.034 1.46E-10 1 0.485 9.45E-06 48 0.182 8.51E-06 
IP5 3 0.037 1E-10 0 0.483 9.23E-06 47 0.232 9.01E-06 
IP6 3 0.021 — 2.03E-10 2 0.500 =: 9.2E-06 9 0.182 7.63E-06 
IP7 3 0.034: 1.8E-10 12 0.490 842E06 48 0.206 9.97E-06 
50,000 IP 3 0.103 4.96E-13 q 1.948 888E-06 50 0.752  8.92E-06 
IP23 0.087 —5.96E-13 5 1.985  7.66E-06 5 0.778 7.09E-06 
IP3 2 0.072 8.87E-06 2 1.910 8.44E-06 48 0.718 1E-05 
IP 3 0.113 4,.96E-13 4 1.949 9.28E-06 50 0.748 — 9.32E-06 
IP5 83 0.101  2.99E-13 3 1.934 9.05E-06 49 0.903  9.87E-06 
IP6—3 0.093 6.95E-13 5 1.978  9.03E-06 5 0.752 — 8.36E-06 
IP7—3 0.079 5.96E-13 5 1.933 8.27E-06 5 0.760 7.65E-06 
100,000. ~IP 2 0.082 6.57E-06 5 4097 9546-06 5 1465 8.83E-06 
P22 0.116 —7.75E-06 6 4347 8.23E-06 52 1.534 —7.02E-06 
IP3 2 0.104 3.14E-06 3 3.826 = 9.07E-06 49 1.392 9.9E-06 
IP 2 0.108  6.95E-06 5 3.872 9.97E-06 5 1.497 9.23E-06 
IP5 2 0.104  4.76E-06 4 3.852 9.73E-06 50 1.473 9.77E-06 
P62 0.095  9.65E-06 6 3.923 9.71E-06 52 1.536 — 8.27E-06 
P72 0.099 8.57E-06 6 3.827 8.89E-06 52 1.537 7.57E-06 


From Tables 3-9, we can observe that the three methods are trying 
to solve (1). However, the improvement and effectiveness of the proposed 
method are pretty straightforward. The tables indicated that modifying 
the IDFDD method in the proposed scheme is a good improvement. The 
MDFDD method remarkably outperforms the IDFDD and IDSL methods 
for nearly all the problems assessed since it has the least number of iter- 
ations, which are far below the number of iterations for the IDFDD and 
IDSL methods. Moreover, the proposed method has less CPU time than the 
IDFDD method. However, the MDFDD method has a higher CPU time than 
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Table 4: Numerical outcomes of MDFDD, IDFDD, and IDSL methods for 
problem 2 


MDFDD IDFDD IDSL 
Dimension IP ITRN  CTM(S) ||/|| ITRN CTM(G) ||K|| ITRN CTM(GS) |All 
100 IP 3 0.026 6.15E-06 26 0.022 7.35E-06 37 0.010  8.96E-06 
IP2 20 0.036 9.95E-06 19 0.014 -9.06E-06 28 0.012 — 9.49E-06 
IP3. 15 0.043 5.55E-06 30 0.022 7.09E-06 43 0.021 7.75E-06 
IP4 8 0.018 4.46E-06 27 0.018 9.54E-06 33 0.017 8.71E-06 
IP5 20 0.024 627E-06 64 0.041 9.19E-06 34 0.009 8.09E-06 
IP6 1 0.017 -5.10E-06 30 0.023 7.61E-06 38 0.012 9.29E-06 
IP7 23 0.043 6.62E-06 39 0.024 6.34E-06 43 0.021  9.82E-06 
1,000. IP 9 0.062 286E-07 28 0.020 9.52E-06 40 0.013 9.72E-06 
IP2 7 0.040 7.93E-06 22 0.025 7.51E-06 32 0.025 7.21E-06 
IP3 20 0.072 5.27E-06 32 0.032 9.18E-06 46 0.028 — 8.4E-06 
IP 6 0.050 6.52E-06 30 0.024 7.91E-06 36 0.027 9.45E-06 
IP5 8 0.062 8.70E-06 71 0.054 7.79E-06 34 0.027 7.27E-06 
IP6 12 0.045 5.73E-06 32 0.020 9.85E-06 42 0.039 7.06E-06 
IP7 24 0.065  6.80E-06 43 0.018  5.72E-06 43 0.032 9.51E-06 
10,000 IPI 33 0.300 5.24E-06 31 0.174 7.89E-06 44 0.153 7.38E-06 
P24 0.345 3.48E-06 24 0.143 9.72E-06 35 0.131 7.82E-06 
IP3 14 0.208 7.43E-06 35 0.129  7.61E-06 49 0.185 9.11 E-06 
IP4 35 0.269 3.35E-06 33 0.174 6.55E-06 40 0.101 7.17E-06 
IP5 7 0.220 6.73E-06 70 0.209 8.45E-06 33 0.119  9.98E-06 
IP6 14 0.207 4.16E-07 35 0.119 817E-06 45 0.141 7.65E-06 
IP7 29 0.350 54E-06 44 0.157 8E-06 46 0.137 7.38E-06 
50,000. IPI 37 0.719 283E-06 33 0.443 7.23E-06 46 0.389 8.08E-06 
IP2 45 0.792 2.38E-07 26 0.339 = 8.91E-06 37 0.334 8.57E-06 
IP3 16 0.594 2.96E-06 37 0.461 6.97E-06 51 0.472 9.99E-06 
IPA 43 1.100 9.23E-07 34 0.433 -9.38E-06 42 0.408 —_7.86E-06 
IP5 0.629 4.38E-06 73 0.825 8.01E-06 33 0.281 — 9.93E-06 
IP6 14 0.515  8.34E-06 37 0.556 7.48E-06 47 0.435. — 8.39E-06 
IP7 32 1.103 8.78E-06 45 0.527  8.09E-06 48 0.447 7.49E-06 
100,000 IPI 36 1114 8695-06 34 0.759 654E-06 47 0.786 8E-06 
IP2 50 L775 8.2E-07 27 0.645 8.06E-06 38 0.650 — 8.48E-06 
IP3 0.595  7.05E-07 37 0.865 9.86E-06 52 0.873  9.89E-06 
IPA 42 1.368  5.53E-06 35 0.760 849E-06 43 0.697 7.78E-06 
IP5 7 1.597  6.23E-06 72 1.480 8.15E-06 33 0.578 — 9.92E-06 
IP6 2 1.756 3.55E-06 38 0.833 6.77E-06 48 0.775 8.3E-06 
IP7 34 2416  6.9E-06 46 1.029 9.52E-06 49 0.823 7.3E-06 
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Table 5: Numerical outcomes of MDFDD, IDFDD, and IDSL methods for 


problem 3 
MDFDD IDFDD IDSL 
Dimension IP ITRN CTM(S) ||Fl| ITRN CTMG)  ||F|| ITRN CTM(S) __|IFrll 

100 IP 6 0.034 852E-06 32 0.012 9.96F-06 37 0.010  8.32E-06 
IP2 12 0.033 9.99E-06 31 0.023 8.84E-06 26 = 0.015.—S«7.7E-06 

IP3 200.044. 2.25E-06 = 350.021 6.42E-06 = 41.—S(0.008 SS 8. ATE-06 

IP4 17  =—0.021._~—«1.76E-06 © 32,—S—«(0.025S «8 45E-06 = 386.017 7.56E-06 

IP5 20 0.045—G.07E-06 34.—Ss(0.012Ss«G.95E-06 = 40.— 0.020 7.81 E-06 

IP6 10 ~—-0.023.—Ss«1.06E-06 -28—Ss«<0.016.—S7.37E-06 © 35—Ss«O.017_—s#8.38E-06 

IP7 17 0.030 7.24E-06 310.023, 7E-06 ~~ 35 0.007__—«8.G1E-06 

1,000 (IPl 22 0.071 6326-06 35 0.047 8.26F-06 40 0.015  9.02E-06 
IP2 15 0.039 4.47E-06 34 0.048 + =7.33E-06 ©=29—S(0.012—S—«8..35E-06 
IP3-26=— 0.074. -2.66E-0737—Ss(0.044._—8.32E-06 © 440.047 9.1 9-06 

IP4 24 0.088 3.26E-06 -35s«.0327.01E-06 = 39—S (0.080 8. 2E-06 

IP5 280.051 7.39E-07 36-—S—«0.049S«9.28E-06 = 43-—(0.040~—«8..83E-06 

IP6 12 0.037 «-3.85E-06 30 ~=—0.026-~—-9.54E-06 = 38-—S—«(02.033«9.09-06 

IP7 13 _-0.082_—4.29E-06_ 33 ———0.036_—7.05E-06 380.028 —_—7.25E-06 

10,000 IPI 26 0.295 6.95F-06 38 0.200 685F-06 43 0. 9.78E-06 
IP2 24 0.345 195E-06 36 0.194 9.5E-06 32 0.168  9.05E-06 
IP337—«0.506. = 34E-06 = 40.181 «6.9E-06 = AT—s«D 9.97E-06 

IP4 27 = 0.319 9.21E-06 37.— «178 s«9.09E-06  42,—s«O.118——8.89E-06 

IP5 -30.—Ss«0.263-«9.76E-06 © 39—Ss«.164.—7.73E-06 © 46—0.166 (9.64 E-06 

IP6 17 ~—-0.331_-—=«9.04E-06 330.136 7.92E-06 4 0.139 9.86E-06 

IP7 210.256 7.66E-06 350.1668. 76E-06 4 0.154 __7.64E-06 

50,000 IPl 29 0.760 6.52E-06 39 0.506 981-06 46 0477 7.5E-06 
IP2 29 0.941 7.04E-06 38 0.520» 8.7E-06 34 0.334 9.92E-06 

IP3 40 = 1.244 5.85E-06 410.529 9.87E-06 50 = 0.618 -—7.65E-06 

IP4 28 0.823 «= 3.9F-06 39 0.540 8.32E-06 4 0.491 9.74B-06 

IPS -35)—s«.871~—3.34E-06 = 41.540 7.08E-06 © 49S (0.575 7.4E-06 

IP6 220.825 G.22E-06 = 35488 7.25E-06 4 0.433 7.56E-06 

IP7 24 «1.087 __—8.32E-06 370.486 7.99E-06 430.446 8.36 E-06 

100,000 IPI] 36 2601 S824E-06 40 0.985 8885-06 47 0.888  7.436-06 
IP2 33 = -2.573.-G.18E-06 39 = «1.1327.88E-06 350.727 ~——«9.82E-06 

IP3 432.423. «9.86E-07 421.382 8.94E-06 5 1.080 7.57E-06 

IP4 34 = 2.529 5.42E-06 © 400.998 7.54E-06 = 450.887 9.65E-06 

IP5 422.712 «9.63E-06 = 42—s«i1.056 = 6. 41E-06 © 500.924. 7.3206 

IP6 201.531 =-2.37E-06 = 36.—Ss(0.904.«G.56E-06 450.862 7.A9E-06 

IP7 24-1634. 5.4E-06 = 380.980 7.23E-06 4 0.825 8.27B-06 
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Table 6: Numerical outcomes of MDFDD, IDFDD, and IDSL methods for 
problem 4 


MDFDD IDFDD IDSL 
Dimension IP ITRN  CTM(S) ||/|| ITRN CTM(G) |K|| ITRN  CTM(GS) |All 

100 IPl 37 0.061 9.93E-06 45 0.018 801E-06 45 0.011  8.46E-06 

IP2 32 0.066 9.00E-06 47 0.029 8AT7E-06 46 0.023 8.83E-06 

IP3 36 0.053 8.51E-06 53 0.032 8.04E-06 54 0.025 9.45E-06 

IP4 37 0.067 7.28E-06 47 0.030 8.34E-06 41 0.019  8.61E-06 

IP5 33 0.056 8.14E-06 49 0.030 9.87E-06 51 0.026 8.19E-06 

IP6 41 0.044 6.93E-06 49 0.022 9.01E-06 52 0.010 — 7.39E-06 

IP7 26 0.029 7.11E-06 41 0.022 8.36E-06 44 0.027 9.03E-06 

1,000. IPI 33 0.084 9.83E-06 46 0.036 TE-05 50 0.040 9.31E-06 

IP2 32 0.090 9.25E-06 48 0.030 9.59E-06 49 0.039 _7.86E-06 

IP3 46 0.102 5.74E-06 53 0.043 7.68E-06 59 0.043 7.05E-06 

IP4 28 0.063 7.2E-06 47 0.052 8.51E-06 47 0.017 — 9.56E-06 

IP5 40 0.080 9.51E-06 52 0.059 9.66E-06 50 0.018 — 8.99E-06 

IP6 37 0.101 849E-06 49 0.032 9.95E-06 55 0.052 7.38E-06 

IP7 24 0.052 7.93E-06 43 0.032 9.41E-06 48 0.038 7.61E-06 

10,000 IPI 36 0.657 9.61E-06 48 0.212 889E-06 53 0.189  9.81E-06 

1P2 40 0.567 3.06E-06 50 0.215 9.42E-06 5 0.227 9.17E-06 

IP3 58 0.594 9.09E-06 53 0.234 = 8.38E-06 64 0.216 — 9.44E-06 

IP4 39 0.761 9.19E-06 50 0.235 -9.73E-06 50 0.190  9.54E-06 

P55 0.718 7.54E-06 54 0.283 1E-05 5 0.244 —-8.49E-06 

IP6 49 0.751 5.9E-06 54 0.254 7.38E-06 58 0.190  7.3E-06 

IP7 3 0.341 -9.06E-06 43 0.221 9.2E-06 48 0.165  8.37E-06 

50,000 IP1 36 2321 659E-06 48 0.698 9.53E-06 56 0.595 7.59E-06 

IP2 39 1581  7.05E-06 53 0.770 9.97E-06 54 0.666 8. 81E-06 

IP3 68 2.331 9.27E-06 55 0.810 8.88E-06 63 0.670 8.29E-06 

IPA 34 1.625  5.01E-06 51 0.752 9.55E-06 5 0.537 8.47E-06 

IP5 53 2.435  9.69E-06 55 0.795 8.83E-06 53 0.638 7.7E-06 

IP6 48 2.340 8.73E-06 54 0.834 9.04E-06 6 0.676 9.28E-06 

IP7 38 1.301 2.77E-06 44 0.655 8.59E-06 50 0.560  8.04E-06 

10,0000 IPr 4 7922 9.92E-06 49 1467 9.825-06. 5: 1.132 9.69E-06 

IP2 46 4.246 9.16E-06 53 1.650 9.8E-06 59 1.311 — 9.21E-06 

IP3 73 4.694 5.2E-06 55 1.908  9.38E-06 63 1.342  9.92E-06 

1P4 46 4.739 6.53E-06 52 1.455 8.6E-06 55 1.113 9.46E-06 

IP5 59 5.184 3.46E-06 56 1.598  8.14E-06 58 1.240  9.03E-06 

IP6 53 3.656 9.1E-06 55 1.586 8.23E-06 6 1.309 — 8.11E-06 

IP7 39 2.125 5.45E-06 44 1418  9.99F-06 5 1.079 7.92E-06 
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Table 7: Numerical outcomes of MDFDD, IDFDD, and IDSL methods for 
problem 5 


MDFDD IDFDD IDSL 
Dimension IP ITRN CTM(S) ||/|| IFRN CTM(G) ||K|| ITRN CTM(G) |All 

100 IP 9 0.026 241E-06 48 0.031  9.76E-06 37 0.023 9.63E-06 
IP2 6 0.015 5.8E-06 45 0.015  8.7E-06 35 0.007 -7.63E-06 

IP3 16 0.037 6.63E-06 53 0.033 8.1E-06 4 0.008 — 7.85E-06 

IP4 3 0.026 127E-06 48 0.036 7.74E-06 37 0.007 7.61E-06 

IP5 4 0.032 297E06 51 0.032 8.6E-06 39 0.012 9.66E-06 

IP6 6 0.007 =7.47E-07 48 0.027 7.65E-06 33 0.014 8.53E-06 

IP7 12 0.029 86E-06 44 0.019 7.78E-06 34 0.012 7.58E-06 

1,000 + IP q 0.042 331E-06 53 0.035. 7.83E-06 4 0.025 7.31E-06 
IP2 0 0.036 5.19E-06 49 0.054 9.18E-06 38 0.019 — 8.27E-06 

IP3. 17 0.041 6.63E-06 57 0.048 8.55E-06 0.039 8.51E-06 

IP 9 0.040 3.13E-06 52 0.060 8.17E-06 40 0.013 8.26E-06 

IP5 9 0.028 4.13E-06 55 0.057 9.45E-06 43 0.036 7.66E-06 

IP6 9 0.039 2.42E-07 52 0.030 8.07E-06 36 0.018 9.25E-06 

x7 i 0.052 945E-06 44 0.040. 7.8E-06 34 0.019 _7.6E-06 

10,000 IP 3 0.237 -5.88E-07 57 0.215 8.26E-06 4 0.108 — 7.93E-06 
IP2 5 0.333 6.48E-06 53 0.171  9.69E-06 0.112 8.97E-06 

IP3 16 0.202 166E-06 61 0.192  9.02E-06 7 0.159 — 9.23E-06 

IP 3 0.150 7.8E-06 56 0.254 8.62E-06 43 0.132 8.95E-06 

IP5 9 0.402 3.93E-06 60 0.225 7.63E-06 46 0.213 8.36E-06 

IP6 13 0.213 3.62E-06 56 0.252 8.52E-06 40 0.093 7.02E-06 

IP7 13 0.211 2.74E-07 44 0.153  7.8E-06 34 0.079 7.6E-06 

50,000 IP 0 0.618 1345-07 60 0.664. 81E-06 46 0.398  8.69E-06 
IP2 2 0.552 4.71E-06 56 0.604 9.51E-06 43 0.349 9.83E-06 

IP3 12 0.465 5.39E-06 64 0.722 8.85E-06 50 0.422 7.08E-06 

IP 7 1.063 2.64E-06 59 0.689  846E-06 45 0.441 — 9.81E-06 

IP5 7 0.961 5.8E-06 62 0.693 9.86E-06 48 0.379 9.17E-06 

IP6 1 0.534 5.8E-06 59 0.655 8.36E-06 42 0.332 7.69E-06 

IP7 12 0.682  1.68E-06 44 0.500  7.8E-06 34 0.287 7.6E-06 

100,000 IP i 1.156 7495-06 61 1.299 8715-06 47 1.125 8.615-06 
IP2 9 0.684 4.11E-06 58 1.145  7.77E-06 44 0.871  9.73E-06 
P3114 1.359 7.52E-06 65 1.336 9.51E-06 51 0.805 7.01E-06 

IP 0 1.360 7E-06 60 1.480  9.09E-06 46 0.873 9.71E-06 

IP5 5 1.204 4.48E-06 64 1.405 8.05E-06 49 0.877 9.08E-06 

IP6 17 1.563 1.77E-06 60 1.206 8.99E-06 43 0.699 7.62E-06 

IP7 1 1.106 1.63E-06 44 1.113 7.8E-06 34 0.520 7.6E-06 
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Table 8: Numerical outcomes of MDFDD, IDFDD, and IDSL methods for 
problem 6 


MDFDD IDFDD IDSL 
Dimension IP ITRN CTM(S) ||/;|| ITRN CTMG) ||K|| ITRN CTMG) ||| 
100 IPl 27 0.205 846E-06 52 0.128 9.42E-06 38 0.073 8.21E-06 
IP2 20 0.195 6.73E-06 47 0.104 7.62E-06 35 0.072 8.98E-06 
IP3) 44 0.264 8.86E-06 60 0.132 9.66E-06 40 0.068 — 8.38E-06 
IPA 31 0.319 7.88E-06 51 0.115  9.37E-06 37 0.068  9.22E-06 
IP5 34 0.264 46E-06 55 0.123 9.69E-06 40 0.071 — 7.5E-06 
IP6 24 0.239 «-5.94E-06 47 0.116 8.84E-06 35 0.061 9.78E-06 
IP7 22 0.213 7.94E-06 45 0.106 7.77E-06 37 0.067 9.35E-06 
1000. +~=IPr 31 2441 9.82E-06 54 154 9.19E-06 41 0.664  8.53E-06 
IP2 25 2.532 8.98E-06 50 092 9.64E-06 38 0.608 — 9.09E-06 
IP3 59 2.687  8.13E-06 57 178 8.11E-06 43 0.697 — 8.2E-06 
IPA 31 2.273 8.28E-06 53 135 8.6E-06 40 0.638 9.51E-06 
IP5 45 2.836 8.66E-06 55 133. 9.24E-06 43 0.680 -8.2E-06 
IP6 25 2.066 8.18E-06 49 047 7.81E-06 38 0.614 — 9.21E-06 
IP7 —_-22 2.083 6.55E-06 45 0.988 7.83E-06 37 0.587 9.37E-06 
2,000 IPl 26 5.618 4725-06 56 3.854 7.96E-06 42 2.135 8.42E-06 
IP2 25 6.920 9.43E-06 52 3.626 8.66E-06 39 2.011  8.95E-06 
IP3 59 8.466 7.28E-06 59 4.029 9.91E-06 44 2.260  8.07E-06 
IP4 27 6.605  8.86E-06 55 3.972 84E-06 41 2.090  9.39E-06 
IP5 43 7.773 7.29E-06 60 223. QAZE-06 44 2.283 8.13E-06 
IP6 27 7.955 6.38E-06 51 3.586 7.86E-06 39 1.978  8.99E-06 


IP7 25 8.424 6.73E-06 45 3.235 7.83E-06 37 1.891 9.37E-06 


Table 9: Numerical outcomes of MDFDD, IDFDD, and IDSL methods for 
problem 7 


MDFDD IDFDD IDSL 
Dimension IP ITRN CTM(S) ||| IFRN CTM(G) ||K|| ITRN CTM(GS) |All 
100 IP 17 0.226 9.74E-06 34 0.080  7.68E-06 3 0.053  9.24E-06 
IP2 23 0.253 5.6E-06 41 0.099 819E-06 36 0.063  9.62E-06 
IP3 (26 0.221 8.98E-06 45 0.102 7.32E-06 40 0.067 7.41E-06 
IP 16 0.175 88E-06 38 0.082 7.76E-06 3 0.063 7.14E-06 
IP5 22 0.174 6.76E-06 42 0.110 9.28E-06 38 0.076 7.42E-06 
IP6 24 0.168  8.55E-06 44 0.100 7.63E-06 39 0.088 8.13E-06 
IP7 21 0.164 5.67E-06 42 0.094 9.17E-06 38 0.073 7.26E-06 
1,000. IP 17 1582 861F-06 35 0.743 831E-06 3 0.497  9.65E-06 
IP2 24 1.553 9.77E-06 44 0.940 9.98E-06 40 0.646 7.28E-06 
IP3 28 1.807 4.01E-06 48 1.012 861E-06 43 0.687 7.64E-06 
IPA 22 1.580 7.03E-06 41 0.870 94E-06 37 0.592  7.62E-06 
IP5 25 1827 7.8E-06 46 1.024 — 8.31E-06 0.657 8.08E-06 
IP6 25 1422 7.42E-06 47 1.046 9.26E-06 42 0.670 8.73E-06 
IP7 28 1849 8.2E-06 46 1.069 851E-06 4 0.655 8.27E-06 
2,000. IP 18 5.262 749E-06 36 2.459 7.66E-06 32 1.656 8.08E-06 
IP2 28 5.929 8.16E-06 46 3.167 _7.43E-06 2121 7.21E-06 
IP3 31 6.115 611E06 49 3.397 8. 82E-06 2.235 7.54E-06 
IP4 (kT 6.253 9.63E-06 42 2.907 9.64E-06 38 1.995  7.54E-06 
IP5 30 7.231 5.2E-06 47 3.298  8.55E-06 42 2.194  8.01E-06 
IP6 31 6.675 4.99E-06 48 3.387 -9.5E-06 3 2.283 8.64E-06 
IP7 25 6.070 5.32E-06 47 3.219 8.77E-06 42 2.188  8.22E-06 
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the IDSL method due to the computation of double direction in the MDFDD 
methods. 

Figures 1—2 display the interpretation of the numerical results of each of 
the three methods using Dolan and Moré [5] performance profiles. We achieve 
this by plotting fraction p(r) of problems for each method within 7 of the 
smallest number of iterations and CPU time. As shown in Figures 1 and 
2, the curves representing the MDFDD method remain above the IDFDD 
and IDSL methods in number iterations. Furthermore, it is above the curve 
representing the IDFDD method for the CPU time. Therefore, the proposed 
method outperforms the IDFDD and IDSL methods in fewer iterations and 
is thus the most efficient method. Finally, from the results in Tables 3-9, it 
is evident that the MDFDD method successfully solves problem (1). 


0.8 MDFDD 


0.7 IDSL 7 


IDFDD 


Figure 1: Performance profile with respect to the number of iterations 


4.2 Application in integral equations 


Chandrasekhar and Breen [3] computed H-equation as the solution of the 
nonlinear integral equation that gives the complete nonlinear equations tech- 
nique. The nonlinear integral equation arising in radiative heat transfer prob- 
lem is given by 

* H(y) 
0 t+Ty 


H(z) =1+ cS H(c) dy, (50) 
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MDFDD 


Figure 2: Performance profile with respect to the CPU time (in second) 


with parameter c € [0,1] and H : [0,1] > R is an unknown function. 


Equation (50) can be written as 


1 
c vH(y) | 
A(#)=]1--= —~dy| =1. dl 
@)=|1-§ [ Bay (51) 
tees ; o fl eH) 3. \* 
By multiplying both sides of (51) with (1 eG lg: ag dy) , we have 
c [' cH(y) = 
F(H)(«) = H(«) (1 | e ay) =0, (52) 
2Jo @+y 


which is called the Chandrasekhar H-equation [23]. However, (52) can be 
discretized by using the midpoint quadrature formula 


| Flw)du = b> S45), (53) 


for pj = (j —0.5)h,0<j<l,andh=+. 


n 


As a result, we have the following system of nonlinear equations: 
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n 
Cc Mit; 
2n jan Mi + pj 


P,(x) = x; 


(54) 


which is known as the discretized Chandrasekhar H-equation that can be 
solved by using some iterative methods. If the initial point zo = (1,1,1,...,1)7, 
then the system in (54) has a solution for all c € (0,1). However, the hardest 
part of the problem (54) is that the Jacobian is singular at c = 1. Therefore 
as c approaches 1, the Jacobian approaches the singularity point. Since our 
method is derivative-free, then it has the advantage to solve problem (54) 
even when c approaches 1. 

To highlight the performance of the MDFDD approach furthermore, we 
conduct some numerical experiments by comparing it with the classical NM, 
IDFDD method [10], and IDSL method [11]. The iteration is also set to 
terminate when ||a,41 — &x|| + ||Fx|| < 107° or when the iterations exceed 
1000, but no point of x, satisfying the stopping criterion. We have tried the 
three methods with the starting point of zp = (1,1,1,...,1)7. Furthermore, 
we use the dimensions (n values) 100 to 20,000 to show the performance of 
each of the three methods. 


Table 10: Numerical results of discretized Chandrasekhar H-equation 


NM IDFDD IDSL MDFDD 

Dimension ITER TIME ITER TIME ITER TIME ITER TIME 

100 2 0.078 - - 38 0.029 13 0.015 

500 5 0.785 - - 41 0.017 14 0.015 

c=0.1 1000 5 3.056 - - 42 0.018 11 0.018 
10000 20 1747 - - 45 0.150 12 0.141 

20000 - - - - 46 0.219 20 0.349 

100 3 0.111 - - 38 0.009 9 0.009 

500 5 0.789 - - 41 0.022 17 0.029 

c=0.9 1000 8 3.729 - - 42 0.026 15 0.024 
10000 - - - - 45 0.231 15 0.215 

20000 - - - - 46 0.216 14 0.377 

100 3 0.111 - - 38 0.012 12 0.011 

500 17 0.874 - - 41 0.032 17 0.020 

c=0.99 1000 8 3.733 - - 42 0.019 12 0.018 
10000 - - - - 45 0.123 11 0.125 

20000 - - - - 46 0.218 13 0.285 

100 3 0.112 - - 38 0.019 13 0.014 

500 7 0.883 - - 41 0.013 16 0.018 

c=0.999 1000 8 3.766 - - 42 0.025 16 0.029 
10000 - - - - 45 0.130 13 0.176 

20000 - - - - 46 0.318 12 0.282 


The numerical results of the methods used to solve Chandrasekhar H- 
equation with different values of parameter c, are shown in Table 10. The 
table clearly indicates that the proposed method outperformed the NM be- 
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cause the NM failed when the number of dimension increased. This is due 
to the fact that as c approaches 1, the Jacobian approaches the singularity 
point. Moreover, the CPU time (in second) of the NM is higher than other 
methods because it solved the Jacobian matrix at each iteration. From Table 
10, we can also observe that the IDFDD method has totally failed because it 
has poor numerical performance, as we have made mentioned earlier in the 
introduction section of this article. Although, the IDSL method solved prob- 
lem (54) completely, but it has more number of iterations than the MDFDD 
method. This shows that our method has effectively solved the discretized 
Chandrasekhar H-equation with the least number of iterations and CPU time. 


5 Conclusion 


In this article, numerical comparisons were made using a set of large-scale 
test problems. Furthermore, Tables 3-9 and Figures 1—2 showed that the pre- 
sented method is practically quite efficient because it has fewer iterations than 
the IDFDD and IDSL methods. Furthermore, we have successfully used the 
proposed method to deal with experiments on the Chandrasekhar H-equation 
in radiative heat transfer. The experiments were carried out and reported 
in Table 10 with different c values, demonstrating a better efficiency for the 
MDFDD method. The numerical results showed that the employed method 
solved the discretized integral equation with fewer iterations and CPU time 
than the NM, IDFDD, and IDSL methods. Future research includes apply- 
ing the MDFDD scheme to solve the discretized three-dimensional nonlinear 
Poisson problem. 
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